Multifractal wavelet filter of natural images 
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Natural images are characterized by the multiscaling 
properties of their contrast gradient, in addition to their 
power spectrum. In this work we show that those properties 
uniquely define an intrinsic wavelet and present a suitable 
technique to obtain it from an ensemble of images. Once this 
wavelet is known, images can be represented as expansions 
in the associated wavelet basis. The resulting code has the 
remarkable properties that it separates independent features 
at different resolution level, reducing the redundancy, and 
remains essentially unchanged under changes in the power 
spectrum. The possible generalization of this representation 
to other systems is discussed. 
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I. INTRODUCTION 

Given the complexity and degree of redundancy of nat- 
ural images, the early visual system had to find good cod- 
ing procedures to represent the visual stimuli internally. 
To achieve this goal, the visual system must have learnt 
the regularities present in the environment where the or- 
ganism lived m. If there are image features that tend 
to appear together, a cell responding quasi-optimally to 
them is rather likely to exist. To find such a representa- 
tion one has first to understand the statistical properties 
of visual scenes common in the environment. In par- 
ticular, the relevance of the second order statistics has 
been pointed out some time ago [0, and internal rep- 
resentations that eliminate these correlations have been 
discussed [pi. However, even if whitening represents an 
improvement of the code, it still leaves much geometrical 
structure that should be dealt with more properly [Q. 
A more systematic study of statistical regularities that 
go beyond the two-point correlations has began rather 
recently ^^. A novel approach to understand the sta- 
tistical properties of natural images has been proposed in 
1^ where the non-gaussian statistics of changes in con- 
trast has been characterized and explained by means of 
a stochastic multiplicative process: Contrast changes at 
a given scale are obtained from those at a coarser scale 
by multiplication with an independent random variable. 
This implies a linear relation between the logarithms of 
the variables at two different scales (this property has 
also been discussed in M, although the existence of a 
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multiplicative process was not noticed). A very rich geo- 
metrical structure has emerged from those studies: con- 
trast changes are organized in such a way that pixels in 
the image can be classified according to the strength of 
the singularities of the contrast gradient [Q . It was also 
checked that the multiplicative process is present in very 
different sets of images and in color natural images [p| . 

In this work we show that when those findings are 
taken into account in an image model there appears a 
very compact representation of the visual world. In par- 
ticular, a wavelet filter that guarantees that variations in 
contrast at different scales are related by a multiplicative 
process can be derived experimentally from a dataset of 
natural images, with the only assumption of the existence 
of such a filter. Furthermore, the code that it defines has 
eliminated a great deal of redundancy. 



II. WAVELET BASIS AND MULTIRESOLUTION 

ANALYSIS 

We start by considering the projection of the contrast 
c{x) on a dyadic wavelet set ^ .j:(a;) = ^(2^a; — k). The 
wavelet Vt -^ focuses on the details of the image at the 
scale r, where r = 2^-', j G Z, at the sampling points 
xq = 2~^k. Here k = (fci,fc2), with ki,k2 G Z ^. 
If ^ defines a wavelet basis, the wavelet projections 
T^c{xq) of the field of luminosity contrast. 



^|c(fo) = r / dx'i'rix- xo)c{x) = r {^r,xo\c) 

(1) 

characterize completely the image. This is the reason 
of the name "multiresolution analysis" : the wavelet pro- 
jection is a description of c{x) at the point xq when the 
image is observed at a variable scale r. This type of anal- 
ysis reaches a compromise between localization in posi- 
tion and in spatial frequency |10| . If the discrete basis is 
complete, then the contrast can be expanded in a wavelet 
basis orthogonal to ^ using the wavelet projections as co- 
efficients. The dual basis ^ ,- verifies that: 

jk 

(*,£l*/fe'> - 2-2^-%,%, (2) 

and the contrast c{x) can then be expressed as: 



c{x) 






(3) 



where a-j: = T^ 'c{2-^k) = 22J(^.-|c). The wavelet 
basis generated by ^ will be called the representation 
basis. Its dual basis, defined by VP will be referred to as 
the analizer basis. Once one of these basis is known, 
the other is completely defined by eq. (0). 



be verified a posteriori, once the wavelet is known. This 
is done in Section IV as follows: first the the analizer 



wavelet ^ is obtained from the representation wavelet 
^. In turn, '^ can be used to evaluate the coefficients 
a r's and from them the rjr's. Once these are known it 
is finally checked that they are indeed scale- independent, 
identically distributed random variables, so checking the 
self-consistency of the multiplicative process model. 



A. Multiscaling in natural images 

In 0] it was experimentally proved that there exists 
a certain class of wavelets, such that some of their low 
order moments are zero, for which: 
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(4) 



(the angular brackets denote the average over an ensem- 
ble of images), a property known as Self-Similarity (SS). 
It was also found Ifl] that the SS exponents Tp have a 
non-trivial dependence on p which can be explained by 
means of a stochastic multiplicative process (see e.g. 
|ll| ) . In terms of the wavelet coefficients a j:, this process 
is expressed as: 
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where [/?] denotes the vector with components given by 
the integer part (rouding down) of those of k. In [||J^] 
eq. (^) was understood in the statistical sense, where 
the variables rj.j^ are statistically independent of the 
a. , ffci ■ The SS exponents Tp define completely the 

distribution of the random variables 77 g and vice- versa. 
Besides, the distribution of the jy's depends only on the 
ratio between the two scales of the wavelet projections 
M. Here this ratio has been fixed at 2 for any pair of 
consecutive scales. This implies that all the ?7 j's are 
identically distributed. 

In this work we go beyond those results by showing 
that natural images posses an intrinsic wavelet for which 
eq. (m is fulfilled point by point. This means that the 
equality holds for any image, at any scale and position 
and that the variables 77 g can be extracted directly from: 
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In order to obtain these variables and to verify their 
statistical properties, one has first to determine this 
wavelet. Under the assumptions that the 77's obtained 
from eq. (0) are scale independent and equally distributed 
variables, the representation wavelet ^P can be experi- 
mentally obtained from a statistical analysis of the im- 
age ensemble. This will be shown in the next section. 
The validity of these two hypothesis on the rj's has to 




FIG. 1. Representation wavelet Vf for the image ensemble. 
The function is represented in gray levels, the darkest points 
indicate the sites where it takes its smallest values. 



III. THE REPRESENTATION MOTHER 
WAVELET 



We suppose that the contrast field obtained from our 
image dataset can be expanded as a superposition of 
wavelets, eq.(||), and that the a g's verify eq. (||), where 
the ?7 jr's are scale independent equally distributed vari- 
ables. Since the images have finite size we will take 
J > where ^Pq^j covers the whole image and it repre- 
sents the mother wavelet, vt'pfj = ^. For the same rea- 
son, the range oi k — (fci,A:2) at the scale j is bounded 
as: ki,k2 = 0,1,. ..,2^-1. 

Averaging c{x) at each point x, it is found that the av- 
erage contrast can be represented as a simple wavelet su- 
perposition: C{x) = (c)(x) oc agg J2j.k\v\ *jfe(^)- Here 
1 77 1 is the first order moment of the distribution of the 
JT^I's, and we have used the assumptions that all the '7j:'s 
have the same marginal distribution and are indepen- 
dent across the scales. By Fourier transforming this field, 
C(/), one easily obtains the Fourier transform of the rep- 
resentation wavelet, ^(/), that reads: 
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FIG. 2. The wavelet 9 along the horizontal (left) and the 
vertical axis (right). 

where A(/) = (1 - e-^'''f'){l - e^^^^). This expres- 
sion is very appeahng. The right hand side compares the 
average contrast at two consecutive scales (related by a 
factor 2), and it expresses that the wavelet is obtained as 
an observation of the scale transformation properties of 
the images. The average \ri\ has an a priori known value: 
\r]\ = i 0. This allows to use eq. (^ with no a priori 
knowledge about the data. 

To obtain experimentally the representation wavelet, 
eq. (m) was applied to a large ensemble of natural im- 
ages. The data were 200 1024 x 1024 images selected at 
random from van Hateren's dataset |12| . Figure n^ shows 
the representation wavelet obtained with this procedure. 
It exhibits two clear features: it is right-left symmetric 
and roughly up-down antisymmetric with a sharp central 
discontinuity (see Figure 0). The first property is ex- 
pected: our world remains statistically unchanged when 
it is reflected from left to right. The rough up-down anti- 
symmetry and the related discontinuity are probably due 
to the sharp contrast between the sky and the ground. 



IV. THE PROJECTION INTO THE 
EXPERIMENTAL BASIS 




FIG. 3. P.d.f.'s of In \rij\ at the scales j — 2 (diamonds), 3 
(crosses), 4 (squares), 5 (x), 6 (triangles) and 7 (crossed x). 

Both correlation coefficients, denoted as pj and p. j- re- 
spectively, take values in [—1, 1] and give a dimensionless 
measure of the degree of statistical dependence. The ob- 
served values of pj are very close to {\pj\ < 10~^, j > 2), 
what confirms rather well the independence of the ?/,^'s 
under changes of scale. On the contrary, the correlation 
coefficients p jare by no means negligible, although they 
are extremely short-ranged: the variables rj .^ and rj .p are 
independent when they are more than one pixel apart; 
after that distance the two-point correlation p .j- decays 
dramatically. It is important to remark that there is no 
need of spatial independence of the ?7j:'s in our wavelet 
model. For this reason, although short-ranged, the ob- 
served dependence is a significant source of information 
about the remaining statistical structure at a fixed reso- 
lution layer. On the other hand, the 77's define a system 
of almost completely independent variables. 

To confirm the validity of the wavelet representation, 
eqs. (|^) and (0), one has still to check that the 77's are 
identically distributed. Figure y exhibits the distribu- 
tion of In 1 77 j I for different scales j. The correspondence 
between them is really very good. 



Once the analizer wavelet is obtained, it can be used 
to evaluate the coefficients a g and from them, using 
eq. (o), the coefficients t] j:. Then, it is possible to check 
whether the rj -^ are scale-independent, identically dis- 
tributed variables. In the affirmative case, this fact self- 
consistently demonstrates that eq. (R) is the intrinsic 
wavelet we are looking for. It is immediate to check that 
the signs of the ry's are independent of their absolute val- 
ues, I77I, and also scale- independent, so it is only necessary 
to verify the independence of the jjyl's. 

The mutual dependence among the |?7 rl's was esti- 
mated by computing the correlation coefficients between 
two \r]\'s. As the number of possible pairs is very large, 
only two types of correlations were considered, which 
should be maximal by construction: those between con- 
secutive scales (j — 1 and j) at equivalent positions, and 
those between consecutive spatial positions (fc and k + d, 
where the components of d take only the values and 1) 
at a fixed scale j. 



V. WAVELET TRANSPARENCY AND 
DECORRELATION 



It has been frequently argued that the the signal that 
arrives to the primary visual cortex has been already 
decorrelated in previous stages of the visual system |g|,|3| . 
In this case, VI would take care of coding more com- 
plex aspects of images. In this regard, the way in which 
the mother wavelet is constructed guarantees a remark- 
able property of transparency to the power spectrum: 
it defines a code that is somewhat independent of the 
second order statistics of the images. The power spec- 
trum of natural images exhibits a power law behaviour, 
S{f) "^ /^^^ '^■' 0- Under the assumption of transla- 
tional invariance, the application of the dccorrclating fil- 
ter to the contrast is equivalent to multiplication in the 
Fourier domain by f^~^ ■ Denoting the decorrelated con- 
trast as Tic(x), it is immediate to see that it has a rep- 
resentation similar to eq. (^: 



Dc(f) 
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where D^ indicates the appUcation of the decorrelat- 
ing operator to the representation wavelet and a'- — 

2^(^~4-'a .jj. Defining now 77'- ~ i^^^-q.^ it is concluded 
that the decorrelated images also posses a random mul- 
tiplicative process. The new representation wavelet D^ 
can be obtained from an ensemble of decorrelated images 
by means of eq. (0) . It is known that the value of the ex- 
ponent e has fluctuations from image to image |lj] . The 
relevance of the wavelet transparency is that if the visual 
stimulus has been decorrelated before it reaches VI, this 
area can use the type of filters described in this paper re- 
gardless of the precise exponent of the power spectrum. 



coefficients 77 jj. If both properties hold, the compact 
code so obtained would be a valuable tool, and the in- 
terpretation of the possible remaining structure at each 
scale very meaningful. 
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VI. CONCLUSIONS 



We have shown that images can be represented as mul- 
tiresolution objects in terms of an appropriate wavelet 
basis ^, in which each resolution level is an independent 
image. One of the advantages of this representation is 
that it is based on the observed properties of the con- 
trast gradient [pP], what in turn leads to an automatic 
reduction of the redundancy. At the same time the spa- 
tial correlations at a given scale are short-ranged, but 
still informative. 

Once the wavelet is known, it is possible to devise a 
compression technique based on the properties of this 
representation. This would have two important advan- 
tages with respect to other wavelet representations (see 
e.g. p5|): first that the scale layers are independent, and 
second that there exists a simple model for the distribu- 
tion of coefi-ccients [g|-§| that can be used in the coding. 

These results have been obtained with the simplest 
wavelet expansion, but the tools presented in this work 
can be taken as the starting point to look for more real- 
istic visual filters. This search should be directed by the 
experimental observation that cells in VI are edge de- 
tectors. From this perspective, the expansion in eq. (0) 
should be generalized to include an orientational degree 
of freedom. Filters of this type have been proposed in 
|l6| and found from an independent component analysis 
of natural images jl^ . These studies should be com- 
bined with the use of overcomplete basis, this introduces 
redundancy but the representation becomes stable under 
small changes in the images |18] 

This analysis could also be carried out over any phys- 
ical system with multiscaling properties, eq. (0) (for 
instance, turbulent flows in fully developed turbulence 
[pT[). For all these systems eq. (M) could be applied to 
obtain the associated representation wavelet ^. However, 
it would be necessary to check the statements of scale in- 
dependence and identical distribution for the projection 
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